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The strong coupling limit (jSgauge = 0) of lattice QCD with staggered fermions enjoys the same 
non-perturbative properties as continuum QCD, namely confinement and chiral symmetry break- 
ing. In contrast to the situation at weak coupling, the sign problem which appears at finite density 
can be brought under control for a determination of the full (/X, T) phase diagram by Monte Carlo 
simulations. Further difficulties with efficiency and ergodicity of the simulations, especially at 
the strongly first-order, low-T, finite-/! transition, are addressed respectively with a worm algo- 
rithm and multicanonical sampling. Our simulations reveal sizeable corrections to the old results 
of Karsch and Mutter Comparison with analytic mean-field determinations of the phase dia- 
gram shows discrepancies of ff{lQ) in the location of the QCD critical point. 
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Figure 1: Phase diagram of strong coupling QCD obtained analytically in the mean-field approximation, 
from [||]. For vanishing quark mass, the first-order transition (solid line) at low temperature turns second- 
order (dashed line) at a tricritical point (TCP). For non-zero quark mass, the first-order transition ends in a 
critical endpoint (CEP), whose trajectory as a function of the quark mass is shown by the dotted line. Also 
displayed are points obtained by the present MC-study: Red dots are transition points for m = 0.1, the red 
circle shows the CEP for T = 1 /2: (m w 0.038, ;U « 0.58). 



1. Introduction 

We consider lattice QCD with Nc = 3 colors and 1 species of staggered fermions (i.e. 4 
continuum flavors) at vanishing gauge coupling. The grand canonical partition function is given by 



with 



Z{m,^) 



+ 2mY,XxXx 



(1.1) 



(1.2) 



where rjjcy = (v = 0) and (— 1)^p<v'^p otherwise^ Like continuum QCD, this model shows 
confinement and chiral symmetry breaking. Therefore, it has been the object of numerous analytic 
studies since the earliest days of lattice QCD, focusing on the mass spectrum ||] and the (/X, T) 
phase diagram [Q, ^,|^, using increasingly refined treatments, all based on the mean-field approxi- 
mation. These investigations are continuing to this day [^]. In contrast, very few numerical studies 
have been performed. Karsch and Mutter [|I|] showed how to express Z(m,/x) as a gas of loops, 
the monomer-dimer-polymer (MDP) ensemble, where the sign problem arising at finite chemical 
potential is very much reduced. They used this formulation to locate the T ^0 finite-/i transition, 
which they found to agree with mean-field predictions. Karsch et al. [||] also found the = finite- 
T transition to be consistent with the expected 0(2) universality class in the chiral limit, turning 
into a crossover for finite quark mass. But a more complete determination of the (/i, T) phase dia- 
gram is lacking. Moreover, Azcoiti et al. [^] reported ergodicity problems with the MDP algorithm 

^ In our notation, all quantities m,jX,T,.. are dimensionless, and should be understood as including the appropriate 
power of the lattice spacing a. 
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Figure 2: Sample configuration with (a) baryonic loop (dashed line) and (b) Di — D2 polymer loop. 

of [|I|], casting some doubt on the T « 0, finite pL results. This is particularly interesting because of 
a mismatch between the critical value of the chemical potential at T 0, }JLc{T « 0,m = 0) = 0.66, 
according to both mean-field ^ and Monte Carlo [|l]], and the baryon mass Mb ~ 3 according to 
mean-field. One would expect /i^ ~ Ms/'i, unless the nuclear interaction is strong. A determination 
of Mb was performed in [10], using conventional HMC at jSgauge = 0. The value agrees closely with 
mean-field. Thus, we now want to determine the (/X, T) phase diagram by Monte Carlo simulations, 
paying special attention to the value of }JLc{T ~ 0). 

The phase diagram can be understood from symmetry considerations. In the chiral limit m = 0, 



the action Eq. (1.2) enjoys the staggered U{\)y.U{\) symmetry 



J(l)v+ie{x)(l)A 



-i(t)v+ie{x)(t>A 



£{X) 



(1.3) 



where U{1)a breaks spontaneously at small ^ and T, giving rise to a chiral condensate o = 
— (X^a^XxXD- Different mean-field approximations all lead to a phase diagram which is quali- 
tatively similar to Fig. |l| with some quantitative differences. The symmetry is restored at large T or 
/I, by a phase transition which is first-order at low temperature, turning second-order at a tricritical 
point (TCP), much like what is expected in real QCD with Nf = 2 massless flavors. For non-zero 
quark mass m where U{\)a is broken explicitly, the second-order transition becomes crossover, 
and the TCP becomes a critical endpoint (CEP). The behaviour of the CEP as the quark mass is 
increased is of particular interest, given recent unexpected findings for Nf = 3 and (2 + 1 ) flavors 
on coarse lattices JTsj]. The dotted line in Fig. [I] shows the trajectory of the CEP in the mean-field 
approximation: it agrees with conventional expectations for real QCD. 

2. Theoretical background 



The usual strategy to deal with the Grassmann fields X^X Eq. ( |1.1[ ) is to integrate them out, 
which yields the customary fermion determinant. At strong coupling, the absence of a gauge action 



allows for an alternative strategy: one integrates over the gauge links U first Jllp. This leads to (for 
Nc = 3): 



:(m,At) = j SixS^X< 



with 



' x,x+v 



3 

I 

k=0 



{Nc-k)\ 



(2.1) 



(2.2) 
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Figure 3: Monte Carlo history of baryon density, for local Metropolis (left) and Metropolis+global worm 
updates (right), for runs with m = 0.025,4^ x 2 at jJ-c w 0.565 using similar CPU time. 



The new degrees of freedom are color singlets: monomers = Y^aXaxXax, dimers D^^xy = 
^{MxMyf,k =1,2,3, and baryons and antibaryons Bj, = ^SahcXaxZbxZcx, = \£ahcXcxXbxXax- 
Moreover, the Grassmann integration generates a "close-packing" constraint: exactly Nc quarks 
and Nc antiquarks must be present at each site. This implies that baryon loops Cg, representing 



n(jc,)))eCB^x5j, in Eq. (gjj), are self-avoiding. It also implies, for the other sites x, that % + = 
Nc, where rix is the number of monomers on x and ni,^ the dimer occupation number (bond number) 
of the hnk connected to x (see Fig. Taking this into account we arrive at the final expression 



I n 

{n^,ni,,CB} b 



(2.3) 



Baryon loops Cg come in two orientations (it) and carry weight 



(2.4) 



where £(Cb) is a sign factor depending on the loop geometry, and ^ > is the winding number 
around the time direction of extent L,. There also exist self-avoiding, non-oriented meson loops, 
which consist of alternating Di and D2 dimers (see Fig. ||(right)). Thus, for a given loop geometry 
C, the partition function should sum over 4 types of loops: two baryon loops with weights w(C, ±), 
and two meson loops related by Di D2, with weight +1. Karsch and MUtter had the idea of 
regrouping these 4 contributions in 2 sets, by associating with each meson loop half of iw(C, +) + 
w{C,—)). In this way, only non-oriented polymer loops C enter in the partition function, with 
weight 

(l + £(C)cosh(3^L,jU)) . (2.5) 



The sign problem, which would have been severe in Eq. ([2.3|), is now much milder. In particular, 
for jLi = one recovers non-negative weights. Moreover, it turns out that the sign problem is very 
much reduced in comparison with that present in the determinant approach, when one integrates 
over fermions first. This makes the study of QCD at large jj. and low T possible. 



3. Algorithmic issues 



The MDP system defined by (2.3), (2.5) was sampled in [|l|] using a local Metropolis algorithm 
which operates on pairs of neighbouring sites and tries to replace two monomers (one on each 
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Figure 4: (a) Baryon density «b versus chemical potential, for systems x Lf at m = 0.1, (b) same, for 
masses m 0,0.025,0.05,0.1 and system 10^ x 2. 

site) by a dirtier on the connecting link or vice versa. Since a monomer carries weight ~ m, this 
prescription is not ergodic in the chiral limit or the infinite-mass limit. Simulations in [|l|] were 
performed over a narrow range of masses. Even so, ergodicity problems were later reported by 
Azcoiti et al.[P], which cast some doubt on the results of M\. Therefore we supplement the local 



Metropolis update above by a worm algorithm [|12|], which was first adapted to strongly coupled 
gauge theories in [13]. In Fig. ^ we show the Monte Carlo history of the baryon density in a 4^ x 2 
system at low quark mass m at the critical n, for both algorithms. The computer time spent is 
similar in both cases. In the Metropolis case, changing the baryon density proceeds via changing 
the monomer density, which is very unlikely. In the worm case, a pair of monomers (the head and 
tail of the worm) is created; then the head is propagated in a succession of nearest-neighbor hops 
until it meets again the tail and annihilates with it, yielding a new contribution to Z having the same 
monomer density, but substantially different baryon density. This allows for efficient simulations 
over the complete range of quark masses. In particular, simulating near or at the chiral limit poses 
no special problem. 

Nevertheless, we still have another difficulty: at low temperature, the finite-/i transition be- 
comes strongly first-order, which makes a correct sampling of the low- and high-density phases 
problematic. To address this issue, we employ Wang-Landau sampling [|l^] to extract an estimator 
for the probability 

P{0,^,m)^ £ d{0-Oik))wk (3.1) 

k={ny,nb,C} 

of a suitable observable O such as the energy density e = -^^logZ or baryon density ns = 
3^7^ log Z. The inverse of the resulting histogram is then used as a weight for multicanonical 
simulations. Note that the weight Wk defined in Eqs. (23) and (2^) allows for reweighting in both, 
jj. and m, once the numbers of monomers and loops with winding number i for each configuration 
are known. Hence, the resulting data from multicanonical sampling can be safely reweighted to 
parameter regions where the corresponding histogram is sufficiently flat. 



4. Numerical results 



We first reproduced the n = 0,T ^ results of |1C] for the chiral condensate, meson and 
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Figure 5: Probability distribution P{g) for masses m — 0.025,0.038,0.05 tuned to criticality on systems 
X 2, L = 8,10,16. A satisfactory data collapse is observed for the middle mass, using Ising critical 
exponents. 



baryon masses. This is a non-trivial consistency check between the two strategies of first integrating 
over the gauge links (this work) or the fermions (followed by HMC) | IC]. 

Then, we fix the quark mass to m = 0.1 and perform a comparison with ||l|], who found a 
phase transition at /X w 0.69 on an 8^ x 4 system, i.e. at temperature T = 1/4. Fig. |4(a)| con- 
firms the problems of ergodicity reported by Azcoiti et al. when using the local Metropolis (data 
taken from ||^). Already on a 4^ system, the simulation remains on one of the two metastable 
branches (density near zero (black triangles) or near the saturation value of 1 baryon per site (black 
squares)), and one is unable to determine the critical value {j.^ of the chemical potential. The Wang- 
Landau/multicanonical approach allows an ergodic sampling of the critical region. The results (red 
stars, continuous line) indicate /i^ « 0.64, which is significantly smaller than the value » 0.69 
of [|l|], presumably obtained from the end of the metastability branch. Fig. 4(a) also shows higher 
temperature (L, = 2) results indicating a smoother transition (as expected) and a slight shift of He 
to smaller values (as opposed to the theoretical predictions of Fig. [I]). 

We now take full advantage of the worm algorithm and explore the chiral limit at fixed tem- 
perature r = 1/2 in Fig. 4(b). As the quark mass is reduced, He shifts to smaller values, and the 
transition becomes first-order, in qualitative agreement with expectations. The order of the tran- 
sition can only be ascertained by a finite-size scaling study, which is illustrated in Fig. |5[ For 
increasing mass values, the 3 panels each show the probability distribution of the chiral condensate 
for 3 spatial volumes. In each case, the distribution is reweighted to the pseudo-critical value of 
jj.. The transition is clearly first-order for the smallest mass, and crossover for the largest. For the 

y 

middle mass, an approximate data collapse is obtained by rescaling the condensate by Lzv, using 
7= 1.237 and v = 0.631 characteristic of the 3d Ising universality class. Indeed, for m > the U(l) 
chiral symmetry is broken explicitly, so we expect no special symmetry breaking other than Z(2) 
to happen at the transition. Thus, we determine the critical endpoint for temperature T = 1/2 to 
be approximately (m^. 0.038, /ic ~ 0.58) (see Fig.[T]). This can be compared with the mean-field 
prediction (wic ~ 0.4, /Xf 0.81) (Fig. |I|) [Q]. The discrepancy is an order of magnitude in m^l 
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5. Outlook and conclusions 

We have presented first results of our study of the strong coupling limit of lattice QCD at fi- 
nite /i and T. The comparison with mean-field results and with the original Monte Carlo study 
of Karsch et al. clearly justifies our project. Algorithmic advances largely suppress the ergodicity 
problems of the latter study, and lead to reliable, new estimates for over a range of masses. 
The large discrepancy between exact Monte Carlo and approximate mean-field determinations of 
the CEP at r = 1/2 emphasizes the need for an exact determination of the whole phase diagram. 
Quantitative mean-field results should be considered with caution. 

Our next step includes the determination of the phase diagram in the chiral limit, and the introduc- 
tion of asymmetric couplings to vary the temperature continuously. As simulations in the chiral 
limit do not pose any problem for a wide range of parameters, further topics of interest might 
include a detailed comparison with chiral perturbation theory, and p ^nn decay. 
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